In [1]:
from pystac_client import Client
import matplotlib.pyplot as plt
import matplotlib
import pandas as pd
import numpy as np
import hvplot.xarray
from hvplot import hvPlot
import xarray as xr
import rioxarray as rxr
import stackstac
import h5py 
import os
import datetime
import geopandas as gpd
import folium
import folium.plugins
from folium.plugins import MousePosition
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.img_tiles as cimgt
import requests
import surfaceAreaGrid
from geopy import geocoders
# For reading COG to xarray
import stackstac
# For displaying image in a jupyter notebook
from IPython.display import Image, display
# For reading from S3 bucket
import boto3
pd.options.plotting.backend = 'holoviews'

Intro: Explore GRA2PES Dataset in the STAC Catalog

In [2]:
stac_api_url = "http://earth.gov/ghgcenter/api/stac/"
In [3]:
# Open catalog, print info and available collections
catalog = Client.open(stac_api_url)
# print catalog and collections titles
print('>> Title: ',catalog.title)
print('>> Collections:')
for collection in catalog.get_collections():
    print(collection.id)

>> Title:  US GHG Center STAC API
>> Collections:
ct-ch4-monthgrid-v2023
emit-ch4plume-v1
lpjwsl-wetlandch4-monthgrid-v1
vulcan-ffco2-yeargrid-v4
goes-ch4plume-v1
odiac-ffco2-monthgrid-v2023
tm54dvar-ch4flux-monthgrid-v1
lpjeosim-wetlandch4-monthgrid-v2
gosat-based-ch4budget-yeargrid-v2
lpjeosim-wetlandch4-daygrid-v2
gra2pes-ghg-monthgrid-v1
gosat-based-ch4budget-yeargrid-v1
micasa-carbonflux-daygrid-v1
sedac-popdensity-yeargrid5yr-v4.11
vulcan-ffco2-elc-res-yeargrid-v4
tm54dvar-ch4flux-mask-monthgrid-v1
casagfed-carbonflux-monthgrid-v3
micasa-carbonflux-monthgrid-v1
lpjeosim-wetlandch4-daygrid-v1
lpjeosim-wetlandch4-monthgrid-v1
oco2geos-co2-daygrid-v10r
epa-ch4emission-yeargrid-v2express
lpjwsl-wetlandch4-daygrid-v1
oco2-mip-meanco2budget-yeargrid-v1
oco2-mip-co2budget-yeargrid-v1
eccodarwin-co2flux-monthgrid-v5
odiac-ffco2-monthgrid-v2022
In [4]:
# Select GRA2PES collection: browse item_assets, spatial and temporal extent
collection = catalog.get_collection('gra2pes-ghg-monthgrid-v1')
collection
Out[4]:
  • type "Collection"
  • id "gra2pes-ghg-monthgrid-v1"
  • stac_version "1.0.0"
  • description "The Greenhouse gas And Air Pollutants Emissions System (GRA²PES) dataset at the GHG Center is an aggregated, regridded, monthly high-resolution (0.036 x 0.036°) data product with emissions of both greenhouse gases and air pollutants developed in a consistent framework. The dataset contains emissions over the contiguous United States covering major anthropogenic sectors, including energy, industrial fuel combustion and processes, commercial and residential combustion, oil and gas production, on-road and off-road transportation, etc. Carbon dioxide (CO₂) emissions are developed along with those of air pollutants including CO, NOₓ, SO₂, and PM2.5 with consistency in spatial and temporal distributions. Emissions by sectors are reported as column totals in units of metric tons per km² per month. Spatial-temporal surrogates are developed to distribute CO₂ emissions to grid cells to keep consistency between greenhouse gases and air quality species. The current version of GRA²PES is for 2021. Long-term emissions and more greenhouse gas species (e.g., methane) are under development and will be added in the future. Description of data processing: Derived from wrfchem_d01 files. Totaled in Z space; native 12-hour (00z and 12z) files averaged to get mean daily emissions; weighted average calculated using weekdy, sundy, satdy to produce average emissions rates for given month. Regridded from native LCC projection and 4000m resolution to EPSG 4326. The native dataset can be accessed from the NIST data center: https://doi.org/10.18434/mds2-3520"
  • links[] 4 items
    • 0
      • rel "items"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1/items"
      • type "application/geo+json"
    • 1
      • rel "parent"
      • href "https://earth.gov/ghgcenter/api/stac/"
      • type "application/json"
    • 2
      • rel "root"
      • href "http://earth.gov/ghgcenter/api/stac/"
      • type "application/json"
      • title "US GHG Center STAC API"
    • 3
      • rel "self"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1"
      • type "application/json"
  • stac_extensions[] 2 items
    • 0 "https://stac-extensions.github.io/render/v1.0.0/schema.json"
    • 1 "https://stac-extensions.github.io/item-assets/v1.0.0/schema.json"
  • renders
    • co
      • assets[] 1 items
        • 0 "co"
      • rescale[] 1 items
        • 0[] 2 items
          • 0 0
          • 1 2
      • colormap_name "spectral_r"
    • nox
      • assets[] 1 items
        • 0 "nox"
      • rescale[] 1 items
        • 0[] 2 items
          • 0 0
          • 1 2
      • colormap_name "greens"
    • so2
      • assets[] 1 items
        • 0 "so2"
      • rescale[] 1 items
        • 0[] 2 items
          • 0 0
          • 1 0.5
      • colormap_name "blues"
    • pm25
      • assets[] 1 items
        • 0 "pm25"
      • rescale[] 1 items
        • 0[] 2 items
          • 0 0
          • 1 0.2
      • colormap_name "purples"
    • dashboard
      • assets[] 1 items
        • 0 "co2"
      • rescale[] 1 items
        • 0[] 2 items
          • 0 0
          • 1 2000
      • colormap_name "spectral_r"
  • item_assets
    • co
      • type "image/tiff; application=geotiff; profile=cloud-optimized"
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
      • title "CO Emissions"
      • description "Estimated monthly total carbon monoxide (CO) emissions across all sectors."
    • co2
      • type "image/tiff; application=geotiff; profile=cloud-optimized"
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
      • title "CO₂ Emissions"
      • description "Estimated monthly total CO₂ emissions across all sectors."
    • nox
      • type "image/tiff; application=geotiff; profile=cloud-optimized"
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
      • title "NOₓ Emissions"
      • description "Estimated monthly total nitrogen oxide (NOₓ) emissions across all sectors."
    • so2
      • type "image/tiff; application=geotiff; profile=cloud-optimized"
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
      • title "SO₂ Emissions"
      • description "Estimated monthly total sulfur dioxide (SO₂) emissions across all sectors."
    • pm25
      • type "image/tiff; application=geotiff; profile=cloud-optimized"
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
      • title "Particulate Matter (PM2.5)"
      • description "Estimated monthly total emissions of fine particulate matter (PM2.5) across all sectors."
  • dashboard:is_periodic True
  • dashboard:time_density "month"
  • title "GRA²PES Greenhouse Gas and Air Quality Species v1"
  • extent
    • spatial
      • bbox[] 1 items
        • 0[] 4 items
          • 0 -128.22655
          • 1 47.89015278
          • 2 -65.30824167
          • 3 22.85824167
    • temporal
      • interval[] 1 items
        • 0[] 2 items
          • 0 "2021-01-01T00:00:00Z"
          • 1 "2021-12-31T00:00:00Z"
  • license "CC-BY-NC-4.0"
  • providers[] 1 items
    • 0
      • name "Colin Harkins"
      • roles[] 2 items
        • 0 "producer"
        • 1 "licensor"
      • url "Colin.harkins@noaa.gov"
In [5]:
# Get items from this collection, examine temporal extent of each
items = collection.get_items()
for item in items:
    print(item)

<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
In [7]:
# Explore one item (one month!): Check out properties, assets
item = collection.get_item('gra2pes-ghg-monthgrid-v1-202112')
item
Out[7]:
  • type "Feature"
  • stac_version "1.0.0"
  • id "gra2pes-ghg-monthgrid-v1-202112"
  • properties
    • end_datetime "2021-12-31T00:00:00+00:00"
    • start_datetime "2021-12-01T00:00:00+00:00"
    • datetime None
  • geometry
    • type "Polygon"
    • coordinates[] 1 items
      • 0[] 5 items
        • 0[] 2 items
          • 0 -137.3143
          • 1 18.173376
        • 1[] 2 items
          • 0 -58.58229999999702
          • 1 18.173376
        • 2[] 2 items
          • 0 -58.58229999999702
          • 1 52.229376000001295
        • 3[] 2 items
          • 0 -137.3143
          • 1 52.229376000001295
        • 4[] 2 items
          • 0 -137.3143
          • 1 18.173376
  • links[] 5 items
    • 0
      • rel "self"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1/items/gra2pes-ghg-monthgrid-v1-202112"
      • type "application/json"
    • 1
      • rel "collection"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1"
      • type "application/json"
    • 2
      • rel "parent"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1"
      • type "application/json"
    • 3
      • rel "root"
      • href "https://earth.gov/ghgcenter/api/stac/collections/gra2pes-ghg-monthgrid-v1"
      • type "application/json"
      • title "GRA²PES Greenhouse Gas and Air Quality Species v1"
    • 4
      • rel "preview"
      • href "https://earth.gov/ghgcenter/api/raster/collections/gra2pes-ghg-monthgrid-v1/items/gra2pes-ghg-monthgrid-v1-202112/map?assets=co2&rescale=0%2C2000&colormap_name=spectral_r"
      • type "text/html"
      • title "Map of Item"
  • assets
    • co
      • href "s3://ghgc-data-store/gra2pes-ghg-monthgrid-v1/GRA2PESv1.0_total_CO_202112.tif"
      • type "image/tiff; application=geotiff"
      • title "CO Emissions"
      • description "Estimated monthly total carbon monoxide (CO) emissions across all sectors."
      • proj:bbox[] 4 items
        • 0 -137.3143
        • 1 52.229376000001295
        • 2 -58.58229999999702
        • 3 18.173376
      • proj:epsg 4326
      • proj:wkt2 "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"
      • proj:shape[] 2 items
        • 0 946
        • 1 2187
      • raster:bands[] 1 items
        • 0
          • scale 1.0
          • nodata -9999.0
          • offset 0.0
          • sampling "area"
          • data_type "float32"
          • histogram
            • max 274.5063781738281
            • min 0.0
            • count 11
            • buckets[] 10 items
              • 0 334968
              • 1 20
              • 2 2
              • 3 2
              • 4 0
              • 5 0
              • 6 0
              • 7 1
              • 8 0
              • 9 1
          • statistics
            • mean 0.15977390873045488
            • stddev 0.9283820264267031
            • maximum 274.5063781738281
            • minimum 0.0
            • valid_percent 73.84708309819413
      • proj:geometry
        • type "Polygon"
        • coordinates[] 1 items
          • 0[] 5 items
            • 0[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
            • 1[] 2 items
              • 0 -58.58229999999702
              • 1 52.229376000001295
            • 2[] 2 items
              • 0 -58.58229999999702
              • 1 18.173376
            • 3[] 2 items
              • 0 -137.3143
              • 1 18.173376
            • 4[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
      • proj:projjson
        • id
          • code 4326
          • authority "EPSG"
        • name "WGS 84"
        • type "GeographicCRS"
        • datum
          • name "World Geodetic System 1984"
          • type "GeodeticReferenceFrame"
          • ellipsoid
            • name "WGS 84"
            • semi_major_axis 6378137
            • inverse_flattening 298.257223563
        • $schema "https://proj.org/schemas/v0.7/projjson.schema.json"
        • coordinate_system
          • axis[] 2 items
            • 0
              • name "Geodetic latitude"
              • unit "degree"
              • direction "north"
              • abbreviation "Lat"
            • 1
              • name "Geodetic longitude"
              • unit "degree"
              • direction "east"
              • abbreviation "Lon"
          • subtype "ellipsoidal"
      • proj:transform[] 9 items
        • 0 0.036000000000001364
        • 1 0.0
        • 2 -137.3143
        • 3 0.0
        • 4 0.036000000000001364
        • 5 18.173376
        • 6 0.0
        • 7 0.0
        • 8 1.0
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
    • co2
      • href "s3://ghgc-data-store/gra2pes-ghg-monthgrid-v1/GRA2PESv1.0_total_CO2_202112.tif"
      • type "image/tiff; application=geotiff"
      • title "CO₂ Emissions"
      • description "Estimated monthly total CO₂ emissions across all sectors."
      • proj:bbox[] 4 items
        • 0 -137.3143
        • 1 52.229376000001295
        • 2 -58.58229999999702
        • 3 18.173376
      • proj:epsg 4326
      • proj:wkt2 "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"
      • proj:shape[] 2 items
        • 0 946
        • 1 2187
      • raster:bands[] 1 items
        • 0
          • scale 1.0
          • nodata -9999.0
          • offset 0.0
          • sampling "area"
          • data_type "float32"
          • histogram
            • max 24279.318359375
            • min 0.0
            • count 11
            • buckets[] 10 items
              • 0 334533
              • 1 293
              • 2 94
              • 3 38
              • 4 17
              • 5 8
              • 6 4
              • 7 3
              • 8 2
              • 9 2
          • statistics
            • mean 26.803277073619228
            • stddev 247.21265875450922
            • maximum 24279.318359375
            • minimum 0.0
            • valid_percent 73.84708309819413
      • proj:geometry
        • type "Polygon"
        • coordinates[] 1 items
          • 0[] 5 items
            • 0[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
            • 1[] 2 items
              • 0 -58.58229999999702
              • 1 52.229376000001295
            • 2[] 2 items
              • 0 -58.58229999999702
              • 1 18.173376
            • 3[] 2 items
              • 0 -137.3143
              • 1 18.173376
            • 4[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
      • proj:projjson
        • id
          • code 4326
          • authority "EPSG"
        • name "WGS 84"
        • type "GeographicCRS"
        • datum
          • name "World Geodetic System 1984"
          • type "GeodeticReferenceFrame"
          • ellipsoid
            • name "WGS 84"
            • semi_major_axis 6378137
            • inverse_flattening 298.257223563
        • $schema "https://proj.org/schemas/v0.7/projjson.schema.json"
        • coordinate_system
          • axis[] 2 items
            • 0
              • name "Geodetic latitude"
              • unit "degree"
              • direction "north"
              • abbreviation "Lat"
            • 1
              • name "Geodetic longitude"
              • unit "degree"
              • direction "east"
              • abbreviation "Lon"
          • subtype "ellipsoidal"
      • proj:transform[] 9 items
        • 0 0.036000000000001364
        • 1 0.0
        • 2 -137.3143
        • 3 0.0
        • 4 0.036000000000001364
        • 5 18.173376
        • 6 0.0
        • 7 0.0
        • 8 1.0
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
    • nox
      • href "s3://ghgc-data-store/gra2pes-ghg-monthgrid-v1/GRA2PESv1.0_total_NOX_202112.tif"
      • type "image/tiff; application=geotiff"
      • title "NOₓ Emissions"
      • description "Estimated monthly total nitrogen oxide (NOₓ) emissions across all sectors."
      • proj:bbox[] 4 items
        • 0 -137.3143
        • 1 52.229376000001295
        • 2 -58.58229999999702
        • 3 18.173376
      • proj:epsg 4326
      • proj:wkt2 "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"
      • proj:shape[] 2 items
        • 0 946
        • 1 2187
      • raster:bands[] 1 items
        • 0
          • scale 1.0
          • nodata -9999.0
          • offset 0.0
          • sampling "area"
          • data_type "float32"
          • histogram
            • max 23.535017013549805
            • min 0.0
            • count 11
            • buckets[] 10 items
              • 0 334546
              • 1 308
              • 2 82
              • 3 27
              • 4 16
              • 5 7
              • 6 3
              • 7 1
              • 8 1
              • 9 3
          • statistics
            • mean 0.041980669686434685
            • stddev 0.23602874947404007
            • maximum 23.535017013549805
            • minimum 0.0
            • valid_percent 73.84708309819413
      • proj:geometry
        • type "Polygon"
        • coordinates[] 1 items
          • 0[] 5 items
            • 0[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
            • 1[] 2 items
              • 0 -58.58229999999702
              • 1 52.229376000001295
            • 2[] 2 items
              • 0 -58.58229999999702
              • 1 18.173376
            • 3[] 2 items
              • 0 -137.3143
              • 1 18.173376
            • 4[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
      • proj:projjson
        • id
          • code 4326
          • authority "EPSG"
        • name "WGS 84"
        • type "GeographicCRS"
        • datum
          • name "World Geodetic System 1984"
          • type "GeodeticReferenceFrame"
          • ellipsoid
            • name "WGS 84"
            • semi_major_axis 6378137
            • inverse_flattening 298.257223563
        • $schema "https://proj.org/schemas/v0.7/projjson.schema.json"
        • coordinate_system
          • axis[] 2 items
            • 0
              • name "Geodetic latitude"
              • unit "degree"
              • direction "north"
              • abbreviation "Lat"
            • 1
              • name "Geodetic longitude"
              • unit "degree"
              • direction "east"
              • abbreviation "Lon"
          • subtype "ellipsoidal"
      • proj:transform[] 9 items
        • 0 0.036000000000001364
        • 1 0.0
        • 2 -137.3143
        • 3 0.0
        • 4 0.036000000000001364
        • 5 18.173376
        • 6 0.0
        • 7 0.0
        • 8 1.0
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
    • so2
      • href "s3://ghgc-data-store/gra2pes-ghg-monthgrid-v1/GRA2PESv1.0_total_SO2_202112.tif"
      • type "image/tiff; application=geotiff"
      • title "Estimated so₂ Emissions"
      • description "Estimated monthly total sulfur dioxide (SO₂) emissions across all sectors."
      • proj:bbox[] 4 items
        • 0 -137.3143
        • 1 52.229376000001295
        • 2 -58.58229999999702
        • 3 18.173376
      • proj:epsg 4326
      • proj:wkt2 "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"
      • proj:shape[] 2 items
        • 0 946
        • 1 2187
      • raster:bands[] 1 items
        • 0
          • scale 1.0
          • nodata -9999.0
          • offset 0.0
          • sampling "area"
          • data_type "float32"
          • histogram
            • max 71.14561462402344
            • min 0.0
            • count 11
            • buckets[] 10 items
              • 0 334881
              • 1 73
              • 2 19
              • 3 11
              • 4 4
              • 5 5
              • 6 0
              • 7 0
              • 8 0
              • 9 1
          • statistics
            • mean 0.016292455984286285
            • stddev 0.35639066244144724
            • maximum 71.14561462402344
            • minimum 0.0
            • valid_percent 73.84708309819413
      • proj:geometry
        • type "Polygon"
        • coordinates[] 1 items
          • 0[] 5 items
            • 0[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
            • 1[] 2 items
              • 0 -58.58229999999702
              • 1 52.229376000001295
            • 2[] 2 items
              • 0 -58.58229999999702
              • 1 18.173376
            • 3[] 2 items
              • 0 -137.3143
              • 1 18.173376
            • 4[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
      • proj:projjson
        • id
          • code 4326
          • authority "EPSG"
        • name "WGS 84"
        • type "GeographicCRS"
        • datum
          • name "World Geodetic System 1984"
          • type "GeodeticReferenceFrame"
          • ellipsoid
            • name "WGS 84"
            • semi_major_axis 6378137
            • inverse_flattening 298.257223563
        • $schema "https://proj.org/schemas/v0.7/projjson.schema.json"
        • coordinate_system
          • axis[] 2 items
            • 0
              • name "Geodetic latitude"
              • unit "degree"
              • direction "north"
              • abbreviation "Lat"
            • 1
              • name "Geodetic longitude"
              • unit "degree"
              • direction "east"
              • abbreviation "Lon"
          • subtype "ellipsoidal"
      • proj:transform[] 9 items
        • 0 0.036000000000001364
        • 1 0.0
        • 2 -137.3143
        • 3 0.0
        • 4 0.036000000000001364
        • 5 18.173376
        • 6 0.0
        • 7 0.0
        • 8 1.0
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
    • pm25
      • href "s3://ghgc-data-store/gra2pes-ghg-monthgrid-v1/GRA2PESv1.0_total_PM25-PRI_202112.tif"
      • type "image/tiff; application=geotiff"
      • title "Estimated PM2.5 Emissions"
      • description "Estimated monthly total emissions of fine particulate matter (PM2.5) across all sectors."
      • proj:bbox[] 4 items
        • 0 -137.3143
        • 1 52.229376000001295
        • 2 -58.58229999999702
        • 3 18.173376
      • proj:epsg 4326
      • proj:wkt2 "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"
      • proj:shape[] 2 items
        • 0 946
        • 1 2187
      • raster:bands[] 1 items
        • 0
          • scale 1.0
          • nodata -9999.0
          • offset 0.0
          • sampling "area"
          • data_type "float32"
          • histogram
            • max 10.905590057373047
            • min 0.0
            • count 11
            • buckets[] 10 items
              • 0 334965
              • 1 20
              • 2 4
              • 3 4
              • 4 0
              • 5 0
              • 6 0
              • 7 0
              • 8 0
              • 9 1
          • statistics
            • mean 0.012796561981333995
            • stddev 0.04822644065182158
            • maximum 10.905590057373047
            • minimum 0.0
            • valid_percent 73.84708309819413
      • proj:geometry
        • type "Polygon"
        • coordinates[] 1 items
          • 0[] 5 items
            • 0[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
            • 1[] 2 items
              • 0 -58.58229999999702
              • 1 52.229376000001295
            • 2[] 2 items
              • 0 -58.58229999999702
              • 1 18.173376
            • 3[] 2 items
              • 0 -137.3143
              • 1 18.173376
            • 4[] 2 items
              • 0 -137.3143
              • 1 52.229376000001295
      • proj:projjson
        • id
          • code 4326
          • authority "EPSG"
        • name "WGS 84"
        • type "GeographicCRS"
        • datum
          • name "World Geodetic System 1984"
          • type "GeodeticReferenceFrame"
          • ellipsoid
            • name "WGS 84"
            • semi_major_axis 6378137
            • inverse_flattening 298.257223563
        • $schema "https://proj.org/schemas/v0.7/projjson.schema.json"
        • coordinate_system
          • axis[] 2 items
            • 0
              • name "Geodetic latitude"
              • unit "degree"
              • direction "north"
              • abbreviation "Lat"
            • 1
              • name "Geodetic longitude"
              • unit "degree"
              • direction "east"
              • abbreviation "Lon"
          • subtype "ellipsoidal"
      • proj:transform[] 9 items
        • 0 0.036000000000001364
        • 1 0.0
        • 2 -137.3143
        • 3 0.0
        • 4 0.036000000000001364
        • 5 18.173376
        • 6 0.0
        • 7 0.0
        • 8 1.0
      • roles[] 2 items
        • 0 "data"
        • 1 "layer"
    • rendered_preview
      • href "https://earth.gov/ghgcenter/api/raster/collections/gra2pes-ghg-monthgrid-v1/items/gra2pes-ghg-monthgrid-v1-202112/preview.png?assets=co2&rescale=0%2C2000&colormap_name=spectral_r"
      • type "image/png"
      • title "Rendered preview"
      • rel "preview"
      • roles[] 1 items
        • 0 "overview"
  • bbox[] 4 items
    • 0 -137.3143
    • 1 18.173376
    • 2 -58.58229999999702
    • 3 52.229376000001295
  • stac_extensions[] 2 items
    • 0 "https://stac-extensions.github.io/raster/v1.1.0/schema.json"
    • 1 "https://stac-extensions.github.io/projection/v1.1.0/schema.json"
  • collection "gra2pes-ghg-monthgrid-v1"

Science topic: Compare urban and rural CO2 season cycles

Let's compare an urban area to a rural area - for example, the Chicagoland area vs. central Illinois.

Read data from STAC into xarray Dataset

In [9]:
search = catalog.search(
    collections='gra2pes-ghg-monthgrid-v1',
    datetime=['2021-01-01T00:00:00Z','2021-12-31T00:00:00Z']
)
# Take a look at the items we found
for item in search.item_collection():
    print(item)

<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
In [10]:
# Read data into an xarray Dataset
ds = stackstac.stack(search.item_collection(),epsg=4326).squeeze()

Quick interactive map to explore the data a little bit...

In [11]:
ds.sel(band='co2')[9,:,:].hvplot(x='x',y='y',cmap='Spectral_r',coastline=True,clim=(0,2000))
Out[11]:
In [12]:
# Take a look at the dimensions of the data
print(ds.sel(band='co2').shape)
print(ds.sel(band='co2').dims)

(12, 947, 2188)
('time', 'y', 'x')

Select areas of interest using GEOJSON shapes and rioxarray

Read in AOIs

GEOJSONs generated via geojson.io

In [44]:
#central_IL = gpd.read_file('./data/central_IL.geojson').to_crs(epsg=4326)
central_IL = gpd.read_file('./data/atlanta.geojson').to_crs(epsg=4326)
chicagoland = gpd.read_file('./data/chicagoland.geojson').to_crs(epsg=4326)

Let's view them on a map!

In [47]:
import cartopy.feature as cf
cmap='Spectral_r'    # colormap definition
vmin=0               # minimum value on colorbar
vmax=1000            # maximum value on colorbar
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
p1 = ds.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
# Rural is green
rural_color = '#00A676'
# Urban is pink
urban_color = '#eb7bc0'
central_IL.plot(ax=ax1,edgecolor=rural_color,linewidth=1.5,facecolor='none')
chicagoland.plot(ax=ax1,edgecolor=urban_color, linewidth=1.5,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)
Out[47]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fa753ad6390>
No description has been provided for this image

Clip CO2 field to AOIs using rioxarray

rio.clip has a few options for selection of grid cells relative to the specified geometry. Here we will use the default functionality. See documentation for more info.

In [48]:
rural_clip = ds.rio.clip(geometries=central_IL.geometry.values)
urban_clip = ds.rio.clip(geometries=chicagoland.geometry.values)
# We also need the latitude values reversed to utilize a function correctly in the next step.
rural_clip = rural_clip.reindex(y=rural_clip.y[::-1])
urban_clip = urban_clip.reindex(y=urban_clip.y[::-1])

Here's what our selections look like:

In [50]:
import cartopy.feature as cf
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
p1 = rural_clip.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
p2 = urban_clip.sel(band='co2')[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
# Rural is green
rural_color = '#00A676'
# Urban is pink
urban_color = '#eb7bc0'
central_IL.plot(ax=ax1,edgecolor=rural_color,linewidth=1,facecolor='none')
chicagoland.plot(ax=ax1,edgecolor=urban_color, linewidth=1,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)
Out[50]:
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fa752b54f20>
No description has been provided for this image

Calculate and apply area weights for lat/lon grid cells

Emissions are in tons km-2 month-1, but our grid cells are in lat/lon space - calculate how many km are in each grid cell, so we can calculate tons month-1 for each AOI.

In [51]:
# This function returns m2, and we need km2. Note the /1000./1000.
rural_weights = surfaceAreaGrid.surfaceAreaGridd(lat_centers=rural_clip.y.values, \
    lon_centers=rural_clip.x.values,ret_area=True)/1000./1000.
urban_weights = surfaceAreaGrid.surfaceAreaGridd(lat_centers=urban_clip.y.values,\
    lon_centers=urban_clip.x.values,ret_area=True)/1000./1000.
In [52]:
# Apply these weights - simple as a multiplication
rural_weighted = rural_clip.sel(band='co2')*rural_weights
urban_weighted = urban_clip.sel(band='co2')*urban_weights
In [53]:
# Take the mean of emissions over AOI grid cells
rural_means = rural_weighted.mean(dim=['x','y'])
urban_means = urban_weighted.mean(dim=['x','y'])
In [54]:
# Reformat some times for plotting purposes
times = [sd[0:10] for sd in rural_means.start_datetime.values]
times = [datetime.datetime.strptime(sd,'%Y-%m-%d').strftime('%b') for sd in times]

Plot urban vs. rural CO2 emissions time series for 2021

In [57]:
fig,ax1 = plt.subplots()
# Rural first
ax1.plot(times[::-1],rural_means.values[::-1],marker='o',color=rural_color,label='Chicago',lw=2.5)
#ax1.set_ylabel('Rural',color=rural_color,fontweight='bold')
# Urban next
#ax2 = ax1.twinx()
ax1.plot(times[::-1],urban_means.values[::-1],marker='o',color=urban_color,label='Atlanta',lw=2.5)
#ax2.set_ylabel('Urban',color=urban_color,fontweight='bold')
plt.legend()
plt.title('Mean urban CO$_2$ emissions, 2021')
Out[57]:
Text(0.5, 1.0, 'Mean urban CO$_2$ emissions, 2021')
No description has been provided for this image
In [ ]: